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Abstract — This paper addresses the prohlem of blind demixing 
of instantaneous mixtures in a multiple-input multiple-output 
communication system. The main objective is to present efficient 
blind source separation (BSS) algorithms dedicated to moderate 
or high-order QAM constellations. Four new iterative batch BSS 
algorithms are presented dealing with the multimodulus (MM) 
and alphabet matched (AM) criteria. For the optimization of 
these cost functions, iterative methods of Givens and hyperbolic 
rotations are used. A pre-whitening operation is also utilized to 
reduce the complexity of design problem. It is noticed that the 
designed algorithms using Givens rotations gives satisfactory per¬ 
formance only for large number of samples. However, for small 
number of samples, the algorithms designed by combining both 
Givens and hyperbolic rotations compensate for the ill-whitening 
that occurs in this case and thus Improves the performance. 
Two algorithms dealing with the MM criterion are presented for 
moderate order QAM signals such as 16-QAM. The other two 
dealing with the AM criterion are presented for high-order QAM 
signals. These methods are finally compared with the state of art 
batch BSS algorithms in terms of signal-to-lnterference and noise 
ratio, symbol error rate and convergence rate. Simulation results 
show that the proposed methods outperform the contemporary 
batch BSS algorithms. 

Index Terms — blind source separation, constant modulus al¬ 
gorithm, multimodulus algorithm, constellation matched error, 
alphabet matched algorithm, Givens and hyperbolic rotations 

I. Introduction 

B lind source separation (BSS) is a fundamental signal 
processing technology that has been intensively used 
in many systems including biomedical, audio and industrial 
applications m. In the context of overdetermined multiple- 
input multiple-output (MIMO) systems, BSS aims to find a 
separation matrix using the received signals and a priori infor¬ 
mation about the statistics or the nature of transmitted source 
signals. Usually, in a communication system, the modulation 
technique being used is known a priori. One can utilize such 
information to recover the same properties in the output signal 
and thus estimate the source signal blindly. Various BSS cost 
functions can be found in literature ||T|, ||2| depending upon the 
types of source signals. Among them, the constant modulus 
(CM) criterion for phase/frequency modulated signals such as 
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PSK/FSK and multimodulus (MM) criterion for QAM signals 
have attracted great interest. 

The CM criterion El restricts the squared modulus of the 
output to be a constant, but such algorithms even work for 
non CM signals. They lead to a number of constant modulus 
algorithms (CMA) used for blind equalization H, Q, blind 
beamforming 0, El and BSS El, ||9l. On the other hand, 
the MM criterion ifTOI takes into account the knowledge of 
square QAM constellation. Its respective cost function deals 
with the real and imaginary parts of the signal separately and 
leads to numerous multimodulus algorithms (MMA) used for 
the application of blind equalization iTol and BSS ifTTl-lfOl. 
MMA outperforms the CMA for the case of square QAM, 
which is used in many modern communication systems such 
as LTE ifTrl and WiMAX ffSl . For such advanced systems 
requiring high data rate, high-order modulations having better 
spectral efficiency are used such as 64-QAM. For these high- 
order modulations, MMA leads to considerable amount of 
residual errors and does not ensure low symbol error rate 
(SFR). Thus, in order to improve the performance of BSS 
algorithms for high-order QAM signals, a number of alphabet 
matched (AM) penalty terms EH-ESl were suggested. All 
of these AM cost functions were found to have good local 
convergence properties and therefore require a good initializa¬ 
tion ifTSll . Thus, alphabet matched algorithms (AMA) should 
be used along with either CMA or MMA, as both of them 
have good global convergence properties. 

Out of numerous CMA solutions, the algebraic one named 
Analytical Constant Modulus Algorithm (ACMA) El provides 
an exact separation in the noise-free case. To overcome the 
drawback of numerical complexity of ACMA, two batch BSS 
algorithms namely Givens CMA (G-CMA) and Hyperbolic 
G-CMA (HG-CMA) were presented in 191, which outperform 
ACMA. The adaptive versions of ACMA and G-CMA were 
presented in lT9l and ll20ll . respectively. Similarly, for the MM 
criterion, an adaptive MMA algorithm was presented in ifTTIl . 
which outperforms the Multi-User Kurtosis (MUK) algorithm 
Ea. Seeing the popularity of ACMA, the same analytical 
approach was used for MM signals and thus an Analytical 
Multimodulus Algorithm (AMMA) was presented in m. 
In terms of blind equalization considering a single source, 
a number of AMA were presented using a combination of 
CM/MM and AM cost functions either in hybrid iflSl or dual 
mode II 22 I . It is shown in 1^ that hybrid and dual mode have 
nearly the same performance. An adaptive blind equalization 


SUBMITTED TO IEEE TRANS. SIGNAL PROCESS., JUNE 2016 


2 


algorithm by combining CM and AM cost functions was 
presented in ll2^ which separates all the sources using multi¬ 
stage cascaded equalizers, where the number of equalizers 
were equal to the number of sources. 

A. Contributions 

In this paper, we propose four new batch BSS algorithms 
utilizing the MM and AM criteria for MIMO systems. The 
major contribution includes the optimization of MM/AM cri¬ 
terion using Givens and hyperbolic rotation parameters for 
the case of multiple sources. Two algorithms are designed by 
minimization of MM and AM cost functions using real Givens 
rotations and named as Givens MMA (G-MMA) and Givens 
AMA (G-AMA), respectively. The other two algorithms are 
designed using both real Givens and hyperbolic rotations and 
thus named as Hyperbolic G-MMA (HG-MMA) and Hyper¬ 
bolic G-AMA (HG-AMA). To the best of our knowledge, this 
is the first paper which presents batch BSS algorithms for MM 
and AM criteria. 

Previously, stochastic gradient techniques were used for 
the minimization of MM and AM criteria o, Qa-ina, 
1241, thus all of these algorithms are adaptive and slow in 
convergence. Therefore, we compare our algorithms with batch 
BSS algorithms designed for CM signals such as ACMA, G- 
CMA and HG-CMA, in terms of signal to interference and 
noise ratio (SINK), symbol error rate (SER), and convergence 
rate. 

B. Paper Organization and Notations 

This paper is organized as follows. In Section |II] the data 
model and a brief overview of BSS principle are presented. 
Section m defines the used criteria as well as Givens and 
hyperbolic rotations. The derivation of proposed algorithms 
G-MMA, HG-MMA, G-AMA and HG-AMA is presented 
in Sections IIVI fVl IVII and IVIII respectively. Section IVIIII 
includes some comments to highlight the important features 
of the proposed algorithms. Simulation results are presented 
in Section |IX] and Section concludes the paper. 

Following are the notations used in this paper, x denotes a 
column vector where its zth entry is denoted by Xi. The real 
and imaginary parts of x are denoted by xr and xi. The matrix 
and its (i, j)-th entry are denoted by X and Xij, respectively. If 
the matrix consists of only real elements then it is represented 
as X. I represents the identity matrix. (.)^ and are used 
to represent matrix/vector transpose and complex conjugate 
transpose, respectively, x denotes the pre-filtered variables, t 
is used to denote s/—l. E\] is the mathematical expectation 
operator and |. | denotes the modulus function. 

II. Problem Formulation 

Consider a MIMO system consisting of Nt sources, each 
having a single antenna element and a receiver equipped 
with an array of antennas. All sources transmit their 
signals over the same band of frequencies. Each transmit¬ 
ted source signal s{i) = SR{i) + LSi{i) is drawn from 
an L-ary square QAM constellation where SR{i),si{i) G 


|±1, ±3,..., ±(v/Z — 1)|- The unknown source signal 

s(i) = [si(/) • • • is passed through a flat fading 

channel represented by an unknown mixing matrix A G 
i^NrxNt Yvjjose elements Omn denotes the channel path be¬ 
tween transmitter n and receiver m. The received signal with 
the added noise can be mathematically represented as 

y(/) = As(/) -I- n(/) (1) 

where n{i) = [ni{i) ■■■ is the white noise 

vector of covariance Here, we assume that the mixing 

matrix A is of full column rank which implies that Nj. > Nf. 

The objective is to recover the source signals s{i) without 
prior knowledge of the channel or without the use of training 
sequences (pilots). This is accomplished using BSS which 
relies on the observation vector y(j) only and also uses some 
source’s structural information. In order to recover the source 
signals (up to a permutation and scaling factors ||9l), vve apply 
a {Nt X Nr) separation matrix W according to 

z(/) = Wy(/) = WAs(/) -f Wn(/) = Gs(/) -f n(/) (2) 

where z{i) = [zi{i) ••• ZNt{i)]^ is the estimated source 
signal vector, G = WA is the {Nt x Nt) global system matrix 
and n(j) = Wn(j) is the filtered noise vector. 

In this paper, we consider batch BSS algorithms in which 
Ns samples of the received signal are collected and then a 
separation matrix is applied on the received data packet Y = 
[y(i) ••• yiNs)], so that O and (|2]| can be rewritten as 

Y = AS -f N, Z = WY (3) 

where Z, S and N are defined in a way similar to the definition 
of Y. In what follows, we seek a separation matrix in the form 
W = VB, where B is a {Nt x Nr) pre-whitening matrix that 
can be computed from a covariance matrix as in |l2l, or simply 
a {Nt X Nr) projection matrix onto the signal subspace (since 
pre-whitening is needed only for the G-MMA/G-AMA but not 
for the HG-MMA/HG-AMA methods). Our main contribution 
lies in designing efficient methods for the computation of the 
matrix V in order to minimize cost functions suitable for high- 
order QAM signals. 

HI. Algorithm Design 

The first step of algorithm design is the selection of a 
suitable cost function. Various cost functions can be found 
in the literature depending upon the properties/types of source 
signals. In the considered case of square QAM signals, we 
have selected the following cost functions. 

A. Cost Functions 

For low order square signals we design MMA using MM 
cost function, however for high-order square QAM signals, we 
design AMA using AM cost function as shown next. 
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1) Multimodulus (MM) Cost Function: For multimodulus 
signals e.g., square QAM, one proposes to estimate the matrix 
V by minimizing the MM criterion defined in IfTOll as 


2 ) Hyperbolic Rotations: The non-unitary Hyperbolic rota¬ 
tion is an (m x m) identity matrix, except for the 

four elements 'Hpp,'Hqq,'Hpq and FLqp given by 


Nt 


'l~Lpq 


cosh( 7 ) e''4 sinh( 7 ) 

Jmma(V) = ^E [(4«(*) - RrY + {zlS) - Rif 

(4) 

'Hqp l~Lqq 


6 “'"^ sinh( 7 ) cosh( 7 ) 


where Rn = Rj = E[|sfl(f)|"^]/E[|sfl(f)p] are dispersion 
constants of the real and the imaginary parts, respectively. This 
cost function was designed such that its minimization can be 
interpreted as fitting the signal into a square shaped signal. 
Thus, it contains structural information of QAM signals and 
also has an inherent ability to restore the phase of the signal. 
Moreover, the MM cost function has several advantages over 
the CM one and leads to: i) faster convergence algorithms 
m, ED, ii) carrier phase recovery 1 ^ . iii) less undesirable 
minima 1291 and iv) ease in hardware implementation lf30l . 

2) Alphabet Matched (AM) Cost Function: Out of the 
variety of AM cost functions mi, ini, USD, m, we have 
selected the one presented in as 

Nt 

Jama(V) = ^ E [g{z,,R{i)) + g{z,g{i))\ (5) 


where 7 G [—r,r] ,r > 0 and /3 G [— 7 r/ 2 , 7 r/ 2 ]. Similar to 
Givens rotations, in the real case, /3 = 0. 

C. Motivation for using Real Givens and Hyperbolic Rotations 

For large number of sources Nt, the difficulty to estimate 

V increases. Thus, to simplify the estimation process, similar 
to Jacobi-like algorithms IMl . ISl, we propose to decompose 

V into a product of Nt{Nt — 1) elementary Givens rotations 
as 

v= n n (9) 

eeps 

where Nsweeps denotes the number of iterations. Parameters 
9 and a are computed in order to minimize the MM criterion 
dUi. Consider a unitary transformation Z = Pp.qY, which 
according to 0 only changes the rows ‘j = p’ and ‘j = f 
of Y so that 

for j ¥=P,q, Zpi = cos{e)y. -h e'-" sin(6<)?/ . 


where g(x) is the constellation matched error (CME) term 
defined as 


Zqi = -e Sin(6»)y^. -P cos{9)y^^ 




g(x) = l-sin^-ix^) 


( 6 ) By omitting the constant terms of Z independent of (9, a), 0 
can be re-written as: 


where n G N and 2d is the minimum distance between alpha¬ 
bet points. The CME in 0 satisfies a number of properties 
that shape the high-order square QAM signals including: i) 
it does not favor alphabet members over others, thus it has 
a uniform behavior, ii) it is locally symmetric around each 
alphabet point, and iii) it places the highest penalty at the 
maximum deviation i.e., the midpoint between two alphabet 
points and does not place any penalty for zero errors i.e., at 
the alphabet points. 

The next step is to devise an efficient method for the 
optimization of the previous cost functions. To guarantee a fast 
convergence with relatively easy implementation, we propose 
to decompose the separation matrix V into a product of 
elementary rotations, similar to Jacobi-like algorithms llMl . 
llT5l . used for matrix diagonalization. Hence V is derived 
using a sequence of Givens and hyperbolic rotations, whose 
parameters are computed by minimizing the MM/AM criteria. 


B. Review of Givens and Hyperbolic (Shear) Rotations 

1) Givens Rotations: The unitary Givens rotation 
Qp,q{9,a) is an (m x m) identity matrix except for the four 
entries Qpp, Gqq,Gpq and Gqp given by 


Qpp 

Gpq 


cos(9) 

sm{9) 

Gqp 

Gqq 


—e“''“ sin( 0 ) 

cos{9) 


where 9 G [— 7 r/ 2 , 7 r/ 2 ] and a G [— 7 r/ 2 , 7 r/ 2 ] are angle 
parameters with a = 0 for the real case. 


pq') 




Zpi,R 


+ {zli,i — 


RiY + {Z^q,,l - Rif 


( 11 ) 


where each term z^^ j, z^^ j equals to g] cos(26»)-h 

5 ! cos(20) cos(2a) + pf cos(20) sin(2a) -I- g^ sin(20) cos(a) -|- 
gl sin(26') sin(Q!) -h 5 ® cos( 2 a) -P gj sin( 2 a) -h gf and gl,j = 
1, • • • , 8 are constant terms depending upon the entries of Y. 
As we can see, further analytical simplification and thus the 
solution of (fTTIi is quite complicated. Similar is the case with 
hyperbolic rotations. These difficulties motivated us to come 
up with a different solution explained belowQ. 


IV. Givens MMA (G-MMA) 

Until now, we have been working in the complex domain 
and to deal with the previously mentioned challenges, we will 
now work in the real domain. Hence, matrix Y is converted 
into a real matrix Y containing real and imaginary parts in 
separate rows as defined in (fT^ . Moreover, a special structure 
of matrix V is introduced and maintained while applying 
the rotations. The transformed real received signal and output 
signal can now be written as Y and Z = VY, respectively, 
where 


Y = 


X,_ 


{2NtxN,), V = 


V/ 


-V, 


(2Ntx2Nt) (12) 


^We have presented this work partly in 021 
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Similarly, S and Z are now {2Nt x Ng) real matrices, which 
can be represented in a way similar to the definition of Y 
in (fTSl l. In order to find the required matrix V, considering 
Lemma 1 of 13^ , the following sequence of real Givens 
rotations are used as a counterpart of (|9]) 


v= n n q{())Qp+Nt,q+Nt p,q+Nt (^) 

Nsweepa i<P,q<Nt 
P/9 

Qq,p+Nt (^) n Gp,P+NXe) (13) 

^<P<Nt 


The rotations Qp^q{9) and Qpj^Mt,q+Nt{(^) are applied suc¬ 
cessively using the same angle parameter (6). Similarly, the 
rotations Qp^q+N^(0) and Qq^p+Nt{&) are applied with another 
angle parameter [6). Note that, these rotations are paired in 
this way to preserve the structure of V given in (ITSli ll^ . 
The rotation Qp^p+Nt(9) is applied to deal with the phase 
shift introduced by the diagonal entries of the mixing matrix 
A. The angle parameters {9) , {9) and [9) are computed is 
such a way to minimize the MM criterion (|4]l, using above 
explained iterative method. For that, we express the MM cost 
function in terms of the angle parameter (9). Now, consider 
a unitary transformation Z = Gp^q^, which according to O 
only changes the rows ‘p’ and ‘g’ of Y so that 


Zji 


Xi j ^ F’ 9, X = cos{9)y^. + sm{9)y^^ 
Zqi = - sin{9)y^. + cos{9)y^^ 


(14) 


Similarly, the rotatioijl Gp+Nt,q+Nt with the same angle 
parameter ( 0 ) modifies the rows ‘p + Nt’ and ‘q + Nt’ in 
a similar way as shown in (fT4l i. Now, (|4]i can be rewritten in 
terms of {9) as (omitting the terms of Z that are independent 
of {9) and assuming for simplicity that Rji = Rj = R) 


Ns 


Jmma{9) = 2 [(4 -rY + Xqi - Rf 

i=l 

+ {^p+Nt,i ~ R) + ~ R) 

Using (O and double angle identities we can write 


(15) 




where 


-if +f 

2 \2.pi ■Lq 


/=[cos(2^) sin(20)] 


——t 

^qi ~ 




(16) 


-pi —qi ::L.pitL.qi 

This allows us to express the first two terms in (fTsT l as 


(17) 


TABLE I: Givens MMA (G-MMA) Algorithm 


Initialization: V = l 2 iVt 

1. Pre-whitening: Y = BY 0{NsN^) 

2. Construct real matrix Y using 03 

3. Givens Rotations: {20NsN^) + O (TVs Y£)/S weep 

for n — 1 . do 

for p = 1 : M do 
for q = p : Nt do 
if p = q then 

a) Compute Qp^p+Nt using I2U . <201 and (?) for 9 (GNa) 

b) Y = 6p,p+jVt Y (iNs) 

c) V = 5p_p^_jV(V 

else 

d) Compute & 6p+jVt,ij+JVt using (19), (20) and 0 

for pme (9) ^ {12Nb) 

e) Y = 5p+jVt.9+AtY {8Na) 

f) V = Qp^q QpJgNt,q + Nt^ 

repeat (d to f) for (p, q + Nt) & {q,P + M) using same 
(9) (20Ns) 

end if 
end for 
end for 
end for 

4. Estimate the complex sources from Y using 0 and (H). 


‘q + Nt’ in ([TtI i. Disregarding the constant terms in (fT^ . we 
can express Jmma{9) as a quadratic form 

Ns 

Jmma{9) = ^ [tiC + V = v'^Tv (19) 

1^1 

The solution v° = [p° that minimizes (O is given by 

the unit norm eigenvector of T corresponding to its smallest 
eigenvalue, so using (fTTl) . we can write 

cos(9) = \ ™tl sin(0) = — = (20) 

Using (l20l i. the computation of Qp^g and follows 

directly from ©. Givens rotations Qp^q+Nt (9) and Gq^p+Nt (^) 
are found similarly and applied successively on Y to compute 
the filtered separation matrix V according to (foT l. The Givens 
rotation Gp,p+Nt{y) for ‘p = q’ can be similarly found by 
following the above explained method. By replacing ‘q’ with 
‘p + Nt in (fTTl) and (fTSl l. the cost function © (with the 
constant terms omitted) can be written as 

Ns 

Jmma{9) = X! = v'^Tv (21) 

i=l 

Hence, the solution v° is the least unit norm eigenvector of 
T and Gp,p+Nt{y) is computed using (l20l) and (|7]i. Matrix 
V is initialized as V = l 2 Arj and the overall algorithm is 
summarized in Table |I] 


iXi -Rf + iX ~ R) 2 v'rt AJ V -h 2 


'2 , '2 
y • V - 

—pi —qi 


-R 


(18) 

Similarly, the terms i ^q+Nt i obtained by 

replacing with L corresponding to indices "p + Nt and 


^For simplicity, we keep the notation Y unchanged even though the matrix 
is modified after each rotation. 


V. Hyperbolic G-MMA (HG-MMA) 

For a small number of samples Ng, the pre-whitening 
operation is not effective and thus the transformed mixing 
matrix A may be far from unitary. In this case, the perfor¬ 
mance of G-MMA deteriorates and thus the J-unitary real 
hyperbolic rotations are applied alternatively along with the 
Givens rotations to overcome this limitation. This results in 
an algorithm named Hyperbolic Givens MMA (HG-MMA). 
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So, now the matrix V can be decomposed into a product of 
elementary hyperbolic rotations Hp^q, Givens rotations Qp^q, 
and normalization transformation Afp^q as follows 

V= n n rp,,(0,7)rp+jv„,+jv,(0,7) 

Ns ■weeps 

'^p,q+Nt (^i i')^q,p+Nt ~7) n ^p^p+nM (22) 

l<p<Nt 

where Tp g = Afp^qQp^qHp^q. Similar to the Givens rotations, 
the hyperbolic rotations Hp^q and 'Hpj^Mt,q+Nt applied 
using the same parameter ( 7 ) while Hp qj^jq^ and T-Lqpj^jq^ 
are applied using another same but opposite parameter ( 7 ) and 
(— 7 ), respectively. We will consider dispersion parameters Rr 
and Rj be equal to 1 and use A/’p,q for normalization. Below 
we give a brief of finding the hyperbolic and the normalization 
transformation parameters to minimize the MM criterion (IHi. 


A. Computation of Hyperbolic and Givens rotations 

Let us consider one hyperbolic transformation Z = Hp^q^, 
which modihes Y according to 


(23) 


Zji for j f- P, q, Zpi = cosh(7)y^. + sinh(7)y^^ 

Zqi = sinh( 7 )y^^ + cosh( 7 )y^. 

Now, using hyperbolic double angle identities we obtain 

" + ^ (C - rI) > 4 - I (C - rI) (24) 


=r^ 

pi : 


where 


u=[cosh( 27 ) sinh( 27 )] 


ri= 


h(y . + y .) y y 

2 2-qi' ^pi2-qi 


Similar expressions can be derived for ^ 


and 0 


that are independent of ( 7 ) yields 


Jmma( 7) = u 


Ns 

r,r] + r,rj 

= u^Ru - 2 u^r 


u — 2u 


Ns 


^ r* + r* 


2 = 1 


u = (R + AJ 2 ) 


Using ( |29] |, the constraint equation results in a 4* order 
polynomial equation 


rT(R + AJ2)-^J2(R + AJ2)-^r = 1 


(30) 


Of the four roots of I®, we use the real valu^l of A that 
results in the minimum value of /i(u, A) with a vector u 
satisfying ui > 0. We then solve for u° = [u°,M 2 ]^ from 
and solve for the hyperbolic sine and cosine of ( 7 ) as 


cosh( 7 ) = 


l + < 


arid sinh( 7 ) = 




(31) 


and 


which allows us to construct the hyperbolic rotations 'Rp,q 
Hp+Nt,q+Nt defined in ([Hll. 

For the remaining hyperbolic rotations Hp^q+Nt 
and Hq^p+Nti the optimization problem in (|26| is 
conducted for the other hyperbolic parameter ( 7 ), 

, /,2 N 

?' ' 


where r, 


2 ( 11 .+ y 


2^+Nt,i' 


v ) —y V 

-p+Nt^i tLqiiLp+Nt,i 


1 T 


RpiR-q+Nt,i 


and 


Then, the modified 


optimization problem is minimized using the same method as 
explained above. This provides the solution u° = 
and the hyperbolic angles are obtained using dM) for 
hyperbolic parameter ( 7 ). The computation of the hyperbolic 
rotations 'Hp.q+NtH) 'Hq^p+Nt{~i) follows directly 

from (1311) and ([8]). Note that these rotations are applied using 
same but opposite hyperbolic angle parameter ( 7 ). 

2) Approximate Solution: In this approach, we will con¬ 
sider the linear approximation of hyperbolic sine and cosine 
around zero given by sinh( 27 ) pz 2 sinh( 7 ) and cosh( 27 ) pz 
cosh( 7 ). Now, let us define the elements of symmetric matrix 

as 


"^(25) 

R and vector r used in 

(EH) 

2 

R = 

rii 

ri2 

q+Nt,i' 


r2i 

r22 

terms 





and 


(32) 


Using (l25l l. (l32l i and neglecting the terms independent of ( 7 ), 
the cost function (|2^ can be rewritten as 


(26) Jmma( 7) = cosh( 47 ) 


rii + r 22 


■ sinh(47)ri2 


where f, = [ 2 W.aVmaJ ' 

optimization problem in (l26l l can be solved using either 
Lagrange multiplier method (exact solution) or by taking 
linear approximation of hyperbolic sine and cosine around zero 
(approximate solution). Both methods are discussed below. 

1) Exact Solution: We consider the constrained optimiza¬ 
tion 

min J^(u) = u^Ru — 2r^u s.t. u^J 2 U = 1 (27) 

U 

where J 2 = diag([l —l]) corresponding to cosh^( 27 ) — 
sinh^( 27 ) = 1. The Lagrangian of dZTl l can be written as 

£(u, A) = u^Ru — 2r^u -f A (u^J 2 U — l) (28) 

The solution of this Lagrangian is given by 


— 2cosh(27)ri — 2sinh(27)r2 (33) 

Setting the derivative of ( l33l l w.r.t ( 7 ) to zero and using the 
previous approximation, we obtain 

sinh( 27 ) (rn + r 22 - ri) - cosh( 27 ) (r 2 - ri 2 ) = 0 (34) 

and thus the solution ( 7 ) is 


1 ^ ( r2-ri2 

7 = - arctanh - 

2 V rii-I-r 22 - ri 


(35) 


(29) 


In a similar way, the hyperbolic rotation parameter ( 7 ) can 
be found using appropriate R and r as explained in section 
IV-All The hyperbolic rotations are computed using ([35]) and 
dH) and applied accordingly as explained in section IV-All 
After applying the hyperbolic rotations, Givens rotations are 
applied in a similar way as explained in section |IV] and then 
normalization rotations are applied as explained below. 

^In the case the set of solutions is empty, we set by default A = 0. 
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TABLE II; Hyperbolic Givens MMA (HG-MMA) Algorithm 


Initialization: V = l 2 A^t 


Subspace projection or approximate pre-whitening 

if Nr > Nt 

0{NsNi) 


1. Create real matrix Y using \\2\ 


2. Hyperbolic, Givens Normalization Rotations: 

{AONsNf) + 

0{NsNt) 


for 71 = \ \ N 


for p = 1 : M do 


iov q = p : Nt do 


if p = q then 


a) Apply Givens rotation using (a to c) of Table Q] (lONs) 

else 


b) Compute Rp^q &i 4ip+jVt,q+jVt using 

1311 and for 

(7) 

(12A) 

c) Y = Aff Y 

(SNs) 

d) V = 'Kp.q'Hp+At.q-i-jVtV 


e) Apply Givens rotation using (d to f) of Table [T| (20^3) 

repeat steps (b to e) for (p, q + Nt) & (q,p + Nt) using 

( 0 , 7 ) ( 0 ,- 7 ), respectively 

(40 A) 

end if 


end for 


end for 


f) Compute A/* using 1371 

(6 A A) 

g) Y = Ay 

(2 A A) 

h) V = Av 


end for 



B. Calculating the normalization transformations 

The normalization is applied to compensate for the dis¬ 
persion parameters Rr and Rj. Let’s consider that we have 
transformed only one row ‘p’ of matrix Y, which corresponds 
to the transformation of rows ‘p’ and ‘p + Nf for matrix 
Y . In this case, the normalization transformation A/"p is an 
identity matrix except for the two diagonal elements J\fpp = 
■N'p+Nt,p+Nt = \ the MM cost function (|4|l (with the 
constant terms omitted) becomes 


Ns 


i7mMa(Ap)— 




Taking the derivative of ( |36] | w.r.t (Ap) and setting the result 
to zero gives optimal normalization parameter 


in ini, ED, iii) it deals with the real and imaginary parts 
of the output signal, separately. Thus, it is relatively easier 
to optimize using real Givens and hyperbolic rotations. In 
this section, G-MMA is used as an initialization followed 
by optimization of AM cost function (with n — 1 for CME 
term in (| 6 ll) using real Givens and hyperbolic rotations, which 
results in algorithms G-AMA and HG-AMA. The combination 
of MMA and AMA is not new and recently used by Labed et 
al. IIJTI for the problem of blind equalization. 

After using G-MMA for the initialization, the matrix V is 
updated using the following Givens rotations 


V — Qp,q+Nt{^)G q ,p+Nt mp,q{0) 

pfq 

ep+jv..,+Jv.(0)V"-i (38) 

where n = uq + I, ..., Nsweeps, where Nsweeps is the number 
of iterations of G-AMA until convergence and ng is the 
number of iterations of G-MMA for initializatioifl Let us 
express the AM cost function in terms of the angle parameter 
(9) which is computed such that Jama( 0) is minimized. Using 
similar derivations as before, one can write 

Ns 

Jama=^ [g (Zpi) + g (Zgi) + g (ip+jVt.i) -f g (4+jVt,i)p9) 

i=l 

where the hrst two terms in ( |39] | can be dehned with n = 1 
in (|6]l as 

g (zp.) = 1 - sin^ { + sin{9)y^^) (^) } 

g {zgi) = 1 - sin^ I sin(6»)y^^ -f cos(6»)y^J } 

and the last two terms are obtained by replacing ‘p’ and ‘q’ 
with ‘p + Nf and ‘q + Nf in ( l40l i, respectively. The bounded 
non-linear optimization problem can now be stated as 


Ap — 




^p+Nt,i 


\ Si=l ^p 


V p 


(37) 


y. 


p+Nt,i 


In our simulations, we observed that the normalization 
rotation is not necessary at each step and can be performed 
only once per sweep. In this case, the diagonal entries of 
matrix Af are Afpp = J\fp+N^^p+N^ = Ap given as in dJTl i 
where 1 < p < HG-MMA is presented in Table HJ 


VI. Givens AMA (G-AMA) 

Eor the design of AM algorithms, the AM cost function 
is selected because of the following reasons: i) it satishes 
all three properties presented in Section IIII-A2I which are 
sufficient conditions to shape the cost function for high-order 
square QAM signals, ii) it is the simplest among all other AM 
cost functions and computationally less expensive. Note that, 
the number of computations in this one is independent of the 
number of alphabet points as opposed to AM cost functions 


mm Jama s.t. 0 £ [-7r/4,7r/4] (41) 

The optimization problem in (l4Tl l can be solved either by using 
MATLAB optimization toolbox that can be termed as ‘exact 
solution’ or by using Taylor series approximation of trigono¬ 
metric functions around zero, which will be referred to as ‘ap¬ 
proximate solution’. This approximation can be justihed using 
Eigure [D which plots the values of AMA cost function Jama 
in (|39] | vs. 6 for some random received pre-whitened signal Y 
after 5 sweeps of G-MMA with Nt = 3,Nr = 5, As = 300, 
SNR = 30dB and normalized 64-QAM constellation. Note that 
the optimum 6° is very close to zero. Thus in the following 
section, we show that for a certain range of 9 close to zero, 
the approximation hts very well with the original values of 
the AMA cost function. 

“^As per observations from the rate of convergence for G-MMA, it converges 
in riQ = 5 for the considered cases. However, one can choose the number of 
sweeps as the one corresponding to an almost flat variation of JTmma- 
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A. Exact Solution 

For the exact solution, the objective function in (|39] | 
is passed to the MATLAB optimization toolbojJI ‘fmin- 
searchbnd’ along with 9o — 0.001 as a starting point and 
bounds 9 G [—7r/4,7r/4], in order to find optimum 9° for 
the minimization of (l4T]) . Givens rotation matrices Qp,q{9°) 
and (0°) are then applied to update V according 

to ( l38b . The remaining Givens rotations Gp,q+Nt{(^) 
Qq,p+Nt{9) can be found similarly by replacing subscripts 
accordingly in ( l39b and ( l40b and then computing optimum 
9°. Then, the separation matrix V is updated again according 
to ( l38b . This process is repeated until convergence. 


B. Approximate Solution 

As observed from Figure [T] the optimum 9° is very close to 
zero, thus the Taylor series approximation around zero can be 
applied. Here, we will consider the approximation up to the 
4* order using the following approximate identities 

q 2 n4 

sin(6») f^9 - cos(6») ~ ^ “ y + ^ (42) 

Now, using the approximation in (l42T i to ‘cos(0)’ and ‘sin(0)’ 
in (l40l i and expanding the terms, results in 




4d2 


2d 


where 


= Air^dry . cos 




— . cos 


d 

' Try 


^ 4,4 ryv^ 

+ n y . cos — f- 
-9* \ d 


■pi 


— Q'lt'^dy y sin 

—pi£-qi 


d 


,3, . ryp^ 


= ird^y 


Try, 


y . sin 

2.qi 


■pi 




(44) 


+ Sit^dy .y . cos 

tLpitLqi 




= Tidy . sin 


ny 


■pi 


—pi 


2,2 ( 
n y . cos — j— 
-9* d 


pi , 

cj = tty^^ sin 




, Cq = 1 + cos 




Similarly g {zqi) can be approximated by 


giZqiY 


48d4 


,9»M 


c^e 


12 # 


d^9^A 


1 


1 


4# 


cr0^+^cr0+-cg\45) 


where all the coefficients are obtained by replacing ‘p’ with ‘q’ 
and ‘q’ with ‘p’ in (l44l) . The 3'^‘* term g (ip+^Vt.i) of ( |39] | has 
the same approximation as given in (|43]) . where the coefficients 
are obtained by replacing ‘p’ with ‘p+Nd and ‘g’ with ‘q+Nt’ 
in (l44li . The last term g{zq+Nt,i) of (l3^ is approximated as 



e 

Fig. 1; JTama vs. 9 for random received pre-whitened signal 
after 5 sweeps of G-MMA with Nt = 3,Nr = 5,Ns = 300, 
SNR = 30dB and normalized 64-QAM constellation. 



Fig. 2: Comparison of exact and approximated G-AMA cost 
functions after 5 sweeps of G-MMA with Nt — 3, = 

5, Ns = 300, SNR = 30dB and normalized 64-QAM con¬ 
stellation 


(EH), where the coefficients are obtained by replacing ‘p’ with 
‘q + Nt" and ‘q with ‘p -f Nt in (l44l i. 

Now, using (031 and (03 in cost function (l39b results in the 
4* order polynomial equation 


Jama^ 


48# 


Ci40V 




where the coefficients in (l46l l are given by 


c, + cf + ^ y G {0, 2,4} 

Ns 


2=1 

Ns 


(47) 


+ cf - 


2 = 1 


Taking the gradient of (03 w.r.t. 9 yields 


6^>yAMA(^) 

d9 


12 # 


(C4# 


4# 






(48) 


where the coefficients are the same as defined in ( 03 . Out of 
the three possible roots of ( 03 , the optimum 9° is selected 
which results in minimum value of f7AMA(6*) in (133 . 

To illustrate that the approximation in (03 is good enough, 
we have compared the original cost function and approximated 
one for a certain range of 9 around zero in Figure ID 

Rotations Qp^q+Nt{^) and Qq^p+Nti^) are similarly found 
by replacing subscripts accordingly and computing optimum 
9°. Then, the rotations are applied successively on Y. 


^Our objective here is just to compute an ‘exact’ solution of E), which 
can be obtained by a linesearch algorithm as well. 
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TABLE III; Givens AMA (G-AMA) Algorithm 


Initialization: V = l2Nt 

1. Pre-whitening: Y = BY 

2. Construct real matrix Y using 02) 

3. Givens Rotations: 

for n = X Ng-ijjeeps do 
if n <= no then 

a) Apply G-MMA as given in Table Q] 
else 

for p = 1 : Yt — 1 do 
for q = p 1 : Nt do 

b) Find optimum (6°) using roots of ED which gives 
minimum value of 03 

c) Compute Qp,q &^Qp-\-Nt,^+Nt using JT) for same {9°) 

d) Y_ — Qp,q ^p-\-N± ,q-\-N±^^ 

e) V = Qp,q Qp+Nt,q-\-Nt^ 

repeat (b to e) for (p,q + Nt) (<J,P + Nt) using same 

in 

end for 
end for 
end if 
end for 


In summary, matrix V is initialized as identity matrix, then 
G-MMA is applied for no = 5 followed by the update of 
matrix V according to (l3^ by applying Givens rotations on 
Y using the above method, until convergence. The overall 
algorithm is summarized in Table |III] 

VII. Hyperbolic G-AMA (HG-AMA) 

As stated earlier, for a small number of samples Ng, J- 
unitary real hyperbolic rotations are applied alternatively along 
with the Givens rotations to overcome the limitation of ill- 
whitening. This results in an algorithm named as Hyperbolic 
Givens AMA (HG-AMA). 

For HG-AMA, first of all G-MMA is used for initialization. 
Then, matrix V is updated iteratively until convergence using 
following hyperbolic 'Hp^q and Givens Qp^q rotations 

V"= n Tp,q+Nt{OA)Tq^p+Ntie,-i) 

Tp,q{e, ^)Tp+Nt,q+Nt {0, 7 )V”-' (49) 

where Tp^q = Qp^gHp^q- Let us express the AM cost function 
in terms of parameter ( 7 ) which is computed such that 
,7 ama( 7) is minimized. Now, using similar derivations as 
before, one can write 

Ns 

Jama=Y id (Zpi) + g (Zqi) + g {zp+Nt,i) + 9 (4-bM.i)f50) 

where the first two terms in (fSOl l can be dehned as 

g {zpi) = 1 - sin^ | (cosh( 7 )y^^ -f sinh( 7 )y^^) | 

g (zqi) = 1 - sin^ I (sinh(7)y^. -f cosh(7)y^.) | 

and the last two terms are obtained by replacing ‘p’ and ‘q’ 
with ‘p + Nt’ and ‘q -j- Nt’ in dSTT i. respectively. 

Figure [3] represents JTama in (l50l l vs. ( 7 ) after 5 sweeps 
of G-MMA with Nt = 3, Nr = 5,Ng = 300, SNR = 30dB 
and normalized 64-QAM constellation. It can be noticed that 
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Fig. 3; JTama vs. 7 for random received pre-whitened signal 
after NSweeps = 5 of G-MMA with Nt = 3,Nr = 5,Ns = 
300, SNR = 30dB and normalized 64-QAM constellation. 


optimum ( 7 °) is very close to zero. Thus, we can apply 
Taylor series approximation of hyperbolic and trigonometric 
functions around zero, in order to find the solution of the 
optimization problem in (l50l l. Next, two possible ways are 
detailed to solve this optimization problem. 

A. Exact Solution 

For the exact solution, the objective function dSOl l is passed 
to the toolbox ‘fminsearch’ with 70 = 0.001 as starting point 
which minimizes (EOll and returns the optimum hyperbolic 
rotation parameter ( 7 °). Tip, 9 ( 7 °) and 'Hp+Nt,q+Nt{l°) are 
then computed and applied to update V according to (l49b . 

For rotations 'Hp^q+Nt{l) and (— 7 ), 1®* and 4* 

terms of the objective function in (l50l) are dehned as 

g {zpt)=cos^ {(cosh( 7 )y^^ + sinh( 7 )y^^^^ 

g {zq+Nt,t)=cot? {(sinh( 7 )y^^ + cosh( 7 )y^^^^ J(^)} 

and the 2"^^ and 3’"^^ terms of (ISOl l are obtained by replacing ( 7 ) 
with (— 7 ) and indices ‘p’ and ‘q+Nt’ with ‘q’ and ‘p+Nt’ in 
(I 52 I 1 . respectively. Now, the modihed objective function is used 
to hnd the optimum ( 7 °). Matrices ELp^q+Nt and ELg^p+Nt are 
then computed using the above explained method and applied 
successively on Y. The process is repeated until convergence. 


B. Approximate Solution 

Again we will use here the Taylor series approximation 
of trigonometric angles given in (l42l) and hyperbolic angles 
around zero up to 4* order, which can be written as 

3 2 4 

sinh( 7 ) 7 + cosh( 7 ) ~ 1 + ^ ^ (53) 

6 


Let’s consider the hrst term of ( fSOl l. which is given in (BTI) as 

g {Zpi) = l-sin^ { ( 0036 ( 7 ) 7 ^ -f sinh( 7 )y^^) } (54) 

Now, applying the hyperbolic angle approximation given in 
(I 53 I 1 to ‘cosh( 7 )’ and ‘sinh( 7 )’ in the argument of sine in 
(I 54 I 1 and expanding the terms, we get 


g (zpi) Ri 1 - sin^ 


{(247* + 2#,,7+127,7 

-|-4y . 7 ^ -f y . 7 ' 

—qi ' —pi ‘ 

Finally, the approximation in (l42l) is used leading to 

9{zpi 


m)} 


(55) 


1 

^4^ 


<7^4 
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7 

Fig. 4; Comparison of exact and approximated H-AMA crite¬ 
rion after no = 5 of G-MMA with Nt = 3,Nr = 5, Ns = 300, 
SNR= 30dB and normalized 64-QAM constellation. 


where 


pi 4 ^4 
Q = TT V . COS 




+ 6TT'^dy .n . sin 

tLpitLqi 


— cos 


ny 


■pi 


— cos 

—pi 


'ny 


■pi 


— nd y . sin — -j— 


PI 3 '6 • 

C 3 = 7T y . Sin 




■qi 


, 2 , . 

y,. Sin ¬ 


es?) 


— Zn'^dii .y . cos 

—pz—gz 


'^y„ 


vi 2 /2 

d; = TT y 


Try 


y . cos 

^qi 


■pi 


.. . 

ndy . sin — -j— 


-pi 


pi ' 

c{ = sin 




Cg* = 1 -f COS 


'^y„ 


Similarly, the other terms g {zqi), g (ip+jVt,i) and g {zq+jsf^^i) 
of ( fSOl l can be approximated as ( l56b . where the coefficients 
are obtained by replacing indices accordingly in (|57] |. 

Now, using (l56l l in ( fSOl l leads to 


Jama" 


1 


■C47^4 


48(i4 

where the coefficients in 




are given by 


N. 


Ci = 


z=l 


(cf 


+ cr + 


^+Nt,' 


+cr^‘ 




(59) 


where I G {0,..., 4}. Taking the gradient with respect to ( 7 ) 
of AMA cost function in (fSSl l. we get 


5 Jama( 7 ) 

dj 


^C47^+^C37^-^C.7-^C4 (60) 


Out of the three possible real roots of dhOl l. the optimum ( 7 °) 
is selected such that Jama ( 7 ) is minimum. 

To illustrate that the approximation in (1581) is good enough, 
we have compared the original cost function and approximated 
one for a certain range of ( 7 ) around zero in Figure |4] 

For the remaining matrices Hp^q+Nt ( 7 ) Tiq^p+Nt (“7)> 
the AM cost function can be written as (l50t and all the 
terms have the same approximation as given in (|56] | with the 
replacement of ( 7 ) with ( 7 ) and indices accordingly. Also, 


TABLE IV; Hyperbolic Givens AMA (HG-AMA) Algorithm 


Initialization: V = l2Nt 

Subspace projection or approximate pre-whitening if Nr > Nt 

1. Construct real matrix Y using 01 ) 

2. Hyperbolic & Givens Rotations: 

for 71 — 1 . do 

if n <= no then 

a) Apply G-MMA as given in Table Q] 
else 

for p = 1 : — 1 do 

for q = p 1 : Nt do 

b) Find optimum ( 7 °) using roots of io) which gives 
minimum value of 

c) Compute 'Hp,q k.'Hp+Nt,q+Nt using )D for same ( 7 °) 

d) Y = 'Hp+Nt,q + Nt'N 

e) V = Hp^q "Kp+JVt.q+Nt V 

f) Apply Givens rotations using (b to e) of Table I till 
repeat steps (b to f) for (p, q + Nt) & (q,p + Nt) using 
( 6 °, 7 °) & ( 6 °,— 7 °), respectively 

end for 
end for 
end if 
end for 


for the 2 "“^ and 3'^'* terms, the sign of coefficients ci and C 3 
are opposite to the one shown in (l56T l. Now, using (l56l l the 
optimization problem in (ISOl l can be written as 


Jama" 


48d4 


C 47 




with 


Ci = ^(cf* + c7'^‘’* + cfZG {0,2,4} 

Z^l 

Ns 


Z=1 

Ns 


(62) 


Cl- 0 ?+"^*’* + cf + 7+"^*’*) 


i=r 


The final solution is obtained by zeroing the gradient of dMT) . 
Once we obtain the solution ( 7 °), matrices 7-Lp^q+Nt{j°) and 
-Hq p+ 7 Vt(— 7 °) are computed using ([ 8 ]). The separation matrix 
V is then updated according to ( |49] |. 

In summary, matrix V is initialized as identity matrix then 
after applying 5 sweeps of G-MMA, matrix V is updated 
according to (l49l l by applying Givens and hyperbolic rotations 
successively on Y using the above explained method, until 
convergence. The overall algorithm is summarized in Table 

HYl 


VIII. Practical Considerations 
We provide here some insight into the proposed algorithms. 


A. Numerical Cost 

Taking into account the structure of the rotation matrices, 
the numerical cost of the proposed algorithms are compared 
with other CMA-like BSS algorithms in terms of the number 
of flops per sweep in Table [V] As can be seen from Table |V] 
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TABLE V; Numerical complexity of different BSS algorithms 


BSS Algorithm 

Complexity Order 

HG-AMA 

HmsN'f + 0(NsNt) 

G-AMA 

70NgN'f + 0(NsNt) 

HG-MMA 

AONgN't + 0{NgNt) 

G-MMA 

20NsNf + 0(NsNt) 

HG-CMA 

■iONgN'f + O(NsNt) 

G-CMA 

IbNgN'f -E O(NgNt) 

ACMA 

0{NsNf) 


the proposed algorithms are much cheaper than ACMA and of 
the same cost order as G-CMA and HG-CMA. Moreover, the 
proposed algorithms have very fast convergence (typically less 
than 10 sweeps) as shown next in the simulation experiments. 
Also, HG-AMA is more expensive but has better performance 
than all the other algorithms as can be observed from the 
simulations results. 

B. Adaptive implementation 

The numerical cost of the designed batch algorithms in¬ 
creases linearly with the sample size Ng. Furthermore, in real 
life environments, systems are time varying and hence the 
separation matrix W has to be re-estimated or updated along 
the time axis. Thus, for slowly time varying systems, this 
update can be obtained by using adaptive estimation methods. 
Utilizing a sliding window technique as in ||9l, one can achieve 
such source separation in an adaptive manner with a numerical 
cost proportional to 0{NsN^) where Ng is the window size 
(instead of total sample size Ng). 

C. Complex implementation 

As shown in section IIII-CI the real matrix representation 
has been introduced to overcome the difficulties encountered 
for the optimization of parameters of complex Givens and 
hyperbolic rotations. However, we can observe that the ob¬ 
tained results can be cast into complex matrix forms using the 
following straightforward relations: 

GpA(^)0p+n,,,+nM± ^ ap.,(0,O)Y 

'^p,q{l)'^p+Nt,q+Nti'y)'^ ^ '^P,q{'lj^)'Y- 

gp,,+NAe)Qq,p+NM± ^ ap.,(0,-|)Y 

'H-p.q+Nti'l)'H-q,p+Nti'l)'^ ^ "^p,g(7j 

where all matrices on left side of (l6^ are real and the right 
ones are complex. Somehow, we have replaced the two degrees 
of freedom of Qp^q(9,oi) (resp. 'Hp^q{'y, P)) by the two free 
parameters 6 and 9 (resp. 7 and 7 ). This way we have avoided 
the non-linear optimization discussed in section IIII-CI 

D. Performance 

The main advantage of the proposed algorithms resides 
in their fast convergence in terms of the number of sweeps 
(typically less than 10 sweeps are needed for convergence) 
and also in terms of sample size (typically Ng = O{10Nt) 
is sufficient for the algorithm’s convergence). Comparatively, 
the ACMA method requires Ng = O{10N^) samples for its 
convergence and standard CMA-like methods need even more 
samples to converge to their steady state. 



Fig. 5: Average SINK of exact and approximate solution of 
HG-MMA vs. SNR for Nt = 5, Nr = 7,Ng = 100 and 
NSweeps = 10 considering both 16-QAM and 64-QAM. 


IX. Simulation Results 


In order to evaluate the performance of the proposed algo¬ 
rithms, simulation results are presented in this section. Due 
to the lack of any batch BSS algorithm dealing with the MM 
criterion, we do the comparison with batch BSS algorithms 
dealing with the CM criterion such as ACMA, G-CMA and 
HG-CMA w.r.t. convergence rate, SER and SINR defined by 


1 

SINR = ^ V SINR 

Nf, ^ 


(64) 


i=i 


with 


SINRj 


_ 


(65) 


where SINRj is the signal to interference and noise ratio at 
the jth output with gij = w^a^, where and blj are the 
ith row vector and jth column vector of separation matrix W 
and mixing matrix A, respectively. R„ is the noise covariance 
matrix and Sj is the (1 x Ng) source vector at jth input. 

We consider a MIMO system consisting of 5 transmitters 
and 7 receivers (Aj = 5, Nr = 7) with the data model 
given in Section |II] Every uncoded data symbol transmitted by 
each source is drawn from 16-QAM, 64-QAM and 256-QAM 
constellations. The resulting signals are then passed through 
matrix A, generated randomly at each Monte Carlo run 
with controlled conditioning and with i.i.d complex Gaussian 
variable entries of zero mean and unity variance. The noise 
variance is adjusted according to specified signal to noise ratio 
(SNR). The results are averaged over 1000 Monte Carlo runs. 


A. Experiment 1: Exact vi. Approximate Solution of HG-MMA 

In Figure |5] we compare the exact and approximate solution 
of HG-MMA in terms of SINR vs. SNR for 16-QAM and 
64-QAM constellations. The number of sweeps Nsweeps and 
samples Ng are set equal to 10 and 100, respectively. We notice 
that both the exact and approximate solutions have the same 
performance for the considered constellations. Therefore, in 
the following simulations for the HG-MMA, we will use the 
approximate solution, as it is cheaper and easier to implement. 
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Fig. 6: Average SINR of HG-MMA and G-MMA vs. SNR for 
different Nsweeps considering Nt = 5,Nr = 7, Ng = 150 and 
16-QAM constellation. 


B. Experiment 2: Finding Optimum Number of Sweeps for G- 
MMA and HG-MMA 

In Figure |6] we examine the effect of the number of sweeps 
on the SINR of the G-MMA and HG-MMA for the case of 
Ns = 150 and 16-QAM. We notice that the performance of 
proposed algorithms increases with the number of sweeps and 
remains almost unchanged after 5 sweeps. So, in the following 
simulations we will fix the number of sweeps to 5. 

C. Experiment 3: Exact vi. Approximate Solution of G-AMA 
and HG-AMA 

Now, we compare the performance of exact and approxi¬ 
mate solutions presented for G-AMA and HG-AMA in terms 
of SINR vs. SNR. Figure |7a] and shows the plots for 
Ns = 200, 64-QAM and Ng = 500, 256-QAM constellations, 
respectively. The number of sweeps Nsweeps is fixed at 10, 
where we used 5 sweeps of G-MMA followed by 5 sweeps 
of AM As. From Figure |7] we notice that both the exact and 
approximate solutions have the same performance. Therefore, 
in the following simulations for the G-AMA and HG-AMA, 
we will use the approximate solution, as it is cheaper and 
easier to implement. 

D. Experiment 4: Finding Optimum Number of Sweeps for 
G-AMA and HG-AMA 

In Figure [8] we examine the effect of the number of sweeps 
Nsweeps on the SINR of the G-AMA and HG-AMA for 
the case of Ns = 200 and 64-QAM. We notice that the 
performance of proposed algorithms increases with the number 
of sweeps and remains almost unchanged after 8 sweeps (5 
G-MMA H- 3 AMA sweeps). So, in the following simulations 
we will fix the number of sweeps to 8. 

E. Experiment 5: Comparison of Rate of Convergence 

In Figure |9] we have compared the convergence rate of the 
the proposed and benchmarked algorithms. The SNR is fixed 
at 30 dB and Ns is selected as 200 and 500 for 64-QAM and 
256-QAM, respectively. It can be noticed that G-AMA and 
HG-AMA converge in 8 sweeps, while all other algorithms 




Fig. 7: Average SINR of exact and approximate solution of 
HG-AMA and G-AMA vs. SNR for Nt = 5, Nr = 7, 

NSweeps — 10- 



Fig. 8; Average SINR of HG-AMA and G-AMA vs. SNR for 
different Nsweeps considering Nt = 5,Nr = 7, Ng = 200 and 
64-QAM constellation. 

converge in 5 sweeps. Even though the proposed algorithms G- 
AMA and HG-AMA require 3 extra sweeps, the performance 
is much better than all the other algorithms. Moreover, the 
performance of HG-MMA and G-MMA is better than the HG- 
CMA and G-CMA. 

F. Experiment 6: Effect of the Number of Samples 

Figure llOal and II Obi show the SINR of our proposed 
algorithms vs. the number of samples Ng for 64-QAM and 
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(b) 256-QAM, Ns = 500 

Fig. 9: Average SINR of iterative batch BSS algorithms vs. 
NSweeps for Nt = b, Np = 7 and SNR = 30dB. 



Number of Samples Ns 
(b) 256-QAM 

Fig. 10: Average SINR of batch BSS algorithms vs. Ng for 
Nt = 5,Np = 7, SNR = 30dB and Nsweeps = 8. 


256-QAM constellations, respectively. The SNR and the total 
number of sweeps Nsweeps are fixed at 30 dB, and 8, 
respectively. We notice that as expected, the larger the number 
of samples the better the performance is. However, we observe 
a threshold point after which the gain is not significant as the 
SINR will be essentially limited by the SNR value. It can be 
seen that the performance of AM algorithms is better than MM 
and CM algorithms. Also, HG-AMA takes the lead among all 
other algorithms. 

G. Experiment 7: Comparison based on SER 

Figure lllal and lllbl depict the SER of AM, MM and 
CM algorithms vs. SNR for the case of 64-QAM and 256- 
QAM constellations, respectively. The number of samples 
Ng = 300 and Ng = 900 are considered for the case of 64- 
QAM and 256-QAM, respectively. As noticed previously, the 
performance of the HG-AMA is significantly better than all 
the other algorithms. Similar to other figures, same pattern 
of performance is observed i.e., the HG-AMA takes the lead 
followed by the HG-MMA, G-AMA, G-MMA and then by the 
HG-CMA, G-CMA and ACMA. By observing these figures, 
we can say that HG-AMA is the only algorithm which works 
well for high-order QAM constellations. 

X. Conclusion 

In this paper, fundamental problems with the physical layer 
for MIMO systems are addressed. The targeted problems 


include channel estimation and blind demixing. Mainly, the 
problem focussed here is to design algorithms for high-order 
QAM signals without using pilot symbols. Four new iterative 
batch BSS algorithms are presented; two of them dealing with 
the MM criterion namely G-MMA and HG-MMA and the 
other two dealing with the AM criterion namely G-AMA and 
HG-AMA. The proposed algorithms are designed using a pre¬ 
whitening operation to reduce the complexity of optimization 
problem, followed by a recursive separation method of unitary 
Givens and J-unitary hyperbolic rotations for the minimization 
of MM/AM criteria. Instead of using complex matrices, a real 
transformation is considered where a special structure of the 
separation matrix in the whitened domain is suggested and 
maintained throughout all transformations. 

The proposed algorithms are mainly designed for the blind 
demixing of MIMO systems involving high-order QAM sig¬ 
nals. Simulation results demonstrate their favorable perfor¬ 
mance as compared to the state of the art algorithms dealing 
with the CM criterion such as G-CMA, HG-CMA and ACMA. 
It is noticed that the G-MMA and G-AMA are cheaper 
and more suitable for large number of samples but in the 
case of small number of samples the HG-MMA and HG- 
AMA should be used. Moreover, out of all the currently 
available batch BSS algorithms and the presented ones, the 
alphabet matched algorithm designed by combining Givens 
and hyperbolic rotations (HG-AMA) is the most efficient one 
for high-order QAM signals such as 64-QAM and 256-QAM. 
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Fig. 11: Average SER of AM, MM and CM algorithms vs. 
SNR for Nt = 5, Nr = 7, Ns weeps = 8 and different Ng. 
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